% USRP Digital Down Converter (DDC)  MatLab Simulation & Analysis

% By : Eng. Firas Abbas
% Edited by Lee Kashka to generalize the function
% USRP DDC consists of 4 stages CIC decimation filter followed by optional
% 31 tap large halfband filter and 7 tap small halfband filter.

function usrpddc(Fs_out)

%initialize everything for no Halfband filters
hb=0;
decim=64e6/Fs_out;
R_cic=decim;

%Determine if halfband filters should be used
if(rem(decim,4)==0)
    hb=2;
    R_cic=decim/4;
else
    if(rem(decim,2)==0)
        hb=1;
        R_cic=decim/2;
    end
end

% Input Signal Sampling Frequency in Hz
% This could be even generalized even more to be a flexible clock rate
Fs_in=64e6;   % 64MHz is the default value

% Number of CIC Sections
N_cic=4;

% HBF Decimation Rate for each
R_hbf=2;

%**************************************************************************
%*************************USRP CIC Filter Analysis*************************
%**************************************************************************

% CIC Input Frequency is 64MHz
Fs_cic=Fs_in;

hcic=mfilt.cicdecim(R_cic,1,N_cic);
% Print Filter Information
info(hcic);

% Ploting CIC Filter Frequency Response
%h=fvtool(hcic,'Fs',Fs_cic);

% Do Gain Compensiation
hgain=dfilt.scalar(1/gain(hcic));
hcicnorm=cascade(hgain,hcic);

% Ploting Filter Frequency Response from 0 MHz to Fs_cic/R_cic MHz
h=fvtool(hcicnorm,'Fs',Fs_cic);
axis([0 Fs_cic/(2*1e6) -120 5]);

%**************************************************************************
%************************* USRP HBF Filter Analysis************************
%**************************************************************************

% HBF Input Frequency is Fs_cic/R_cic MHz
Fs_hbf=Fs_cic/R_cic;

%coeff for large hb generated from round(2^18 * halfgen4(.7/4,8))./131072
Numerator = [-0.000816345214843750 0 0.00339508056640625 0 ...
             -0.00969696044921875 0 0.0225753784179688 0 ...
             -0.0465927124023438 0 0.0911941528320313 0 ...
             -0.188491821289063 0 0.628349304199219 1 ...
             0.628349304199219 0 -0.188491821289063 0 0.0911941528320313 ...
             0 -0.0465927124023438 0 0.0225753784179688 0 ...
             -0.00969696044921875 0 0.00339508056640625 ...
             0 -0.000816345214843750];  
         
hbf = dfilt.dffir(Numerator);
info(hbf);
num  = get(hbf, 'Numerator');  % Get the numerator from the current filter.
hhbf=mfilt.firdecim(R_hbf,num);

%coeff for small hb generated from 2 * halfgen4(.75/8,2)
NumSmall = [-0.0815591197584173 0 0.578373182212901 1 0.578373182212901 ...
            0 -0.0815591197584173];
hbfs = dfilt.dffir(NumSmall);
info(hbfs);
num = get(hbfs, 'Numerator');  % Get the numerator from the current filter.
hhbfs=mfilt.firdecim(R_hbf,num);

%Gain Compensiation by dividing by 2
hgain2=dfilt.scalar(1/2);

if(hb==2)
    hhbfnorm=cascade(hgain2,hhbfs,hgain2,hhbf);
    % Plotting Filter Frequency Response
    h=fvtool(hhbfnorm,'Fs',Fs_hbf);
else
    if(hb==1)
       hhbfnorm=cascade(hgain2,hhbf);
       % Ploting Filter Frequency Response
       h=fvtool(hhbfnorm,'Fs',Fs_hbf);
    end
end



%**************************************************************************
%****************************Complete USRP DDC*****************************
%**************************************************************************

if(hb~=0)
    hddc=cascade(hcicnorm,hhbfnorm);
    h=fvtool(hcicnorm,hhbfnorm,hddc,'Fs',[Fs_cic,Fs_hbf,Fs_cic]);
    axis([0 Fs_cic/(1e6*R_cic*hb*2) -98 2]);
    legend(h,'CIC','HBF','CIC+HBF');
else
    hddc=hcicnorm;
    h=fvtool(hddc,'Fs',Fs_cic);
    axis([0 Fs_cic/(1e6*R_cic) -98 2]);
    legend(h,'CIC');
end

end